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Abstract 

The dissipation of kinetic and magnetic energy in the interstellar medium (ISM) can proceed through 
viscous, Ohmic or ambipolar diffusion (AD). It occurs at very small scales compared to the scales at 
which energy is presumed to be injected. This localized heating may impact the ISM evolution but 
also its chemistry, thus providing observable features. Here, we perform 3D spectral simulations of 
decaying magnetohydrodynamic turbulence including the effects of AD. We find that the AD heating 
power spectrum peaks at scales in the inertial range, due to a strong alignment of the magnetic and 
current vectors in the dissipative range. AD affects much greater scales than the AD scale predicted 
by dimensional analysis. We find that energy dissipation is highly concentrated on thin sheets. Its 
probability density function follows a lognormal law with a power-law tail which hints at intermittency, 
a property which we quantify by use of structure function exponents. Finally, we extract structures of 
high dissipation, defined as connected sets of points where the total dissipation is most intense and we 
measure the scaling exponents of their geometric and dynamical characteristics: the inclusion of AD 
favours small sizes in the dissipative range. 


1 Introduction 


In partly ionized astrophysical fluids, magnetic fields remain attached to the charged particles. The 
neutral fluid does not feel the Lorentz force and usually drifts with respect to the charges. Only the 
ion-neutral drag can remind the neutrals about the existence of magnetic fields. The fields are then able 
to slip through the bulk of the fluid: this process is called ambipolar diffusion (AD). 

Mestel and Spitzer |1956| were the first to realize its importance in the context of star formation, 


where it would help magnetic fields to diffuse out of a contracting dense core and allow it to form a star. 
Mullan |197l] discovered how ambipolar diffusion could influence the dynamics of shocks which then 


yielded a wealth of papers on the chemical signatures of C-type shocks, starting with Draine et al. [1983] 


and Flower et al. [1985] . Toth [1995] produced the first multi-dimensional simulations with AD for the 


stability of such shocks, and this opened the way to a collection of analytical and numerical studies in 
various astrophysical contexts. Brandenburg and Zweibel [1994 and Brandenburg and Zweibel [1995] 
envisaged that ambipolar diffusion could form very sharp structures, which would then induce the Ohmic 
resistivity to reconnect the magnetic field in the interstellar medium (ISM): they argued this could be 
a key element in solving observational problems with the galactic dy namo [Zweibel and Brandenburg[ 
1997| . In discs Blaes and Balbus [1994 , Brandenburg et al. [1995 and |Mac Low et al. [1995] showed AT) 
was able to modulate the magneto-rotational instability. In the context of clouds and stars formation 


Nakamura and Li [2005] , Nakamura et al. 2008] , Kudoh and Basu [2008 and Kudoh and Basu [2011] 


computed how the magnetic support of clouds can leak out to let the gas condense and form dense cores 
and stars. Finally, in the context of ISM turbulence Padoan et al. [2000 , Zweibel 2002 , Oishi and Mac 
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Low [2006] , Li et al. [2008] , McKee et al. 2010 , Li et al. [2012 and Li et al. [2012] focused on how 


magnetic fields decouple from the neutrals velocity or density and estimated the heating resulting from 
the ion-neutral drift. 

In the diffuse interstellar medium (ISM), turbulent energy dissipation can be an important source 
of suprathermal energy driving hot chemistry [Falgarone and Puget 1995 . This may be evidenced by 
the observed high values of the column density of species such as CH + and SH + . The formation of 
such species requires energy barriers of the order of at least 2000 K to be overcome, in clouds where the 
average temperature is known to be a few tens of K. One possible explanation is that these cold clouds 
contain pockets of hot gas heated by intermittent turbulent dissipation. Hot chemistry is activated there, 
and it is possib l e to construct models of turbulent dissipation that account for the high column densities 


of CH + Godard et al. 


2009 


In a turbulent magnetized fluid, dissipation can of course be due to viscosity or resistivity, but when 
the fluid is partially ionized and AD is at play, there can also be a significant contribution from the heat 
released by ion-neutral friction, as demonstrated by several authors in the context of the ISM |Scalo[ 


19771 |Zweibel a nd Josafatsson| 


1983 


Elmegreen 


1985 


Padoan et al. 


2000 


Li et al. 


2012 


Not only 

does the heating help to raise the temperature which increases the rate of some chemical reactions, but 
the ion-neutral drift velocity provides additional energy in the reaction frame for ion-neutral reactions. 
In some instances, this can open new chemical routes which would otherwise be blocked by reaction 
barriers. The places of strong AD heating are thus expected to bear specific chemical signatures such as 
the ones encountered in magnetized vortices [Godard et al. 2009] or C-shocks [Lesaffre et al. 2013 . We 
would hence like to characterize the geometry and statistical properties of the regions of strong turbulent 
dissipation in the ISM. 

Before us, Uritsky et al. 20To] conducted a thorough study of the statistics of strong dissipation in 
the context of incompressible pure magnetohydro dynamic (MHD) turbulence. In this paper we make 
some progress towards the physics of the ISM and we work with incompressible MHD turbulence with 
or without AD. We stay within the model of incompressible MHD in a first step to link our work with 
Uritsky et al. [2010] and to allow the use of spectral methods which are well suited for the study of small- 


scale dissipative structures because of their very low level of numerical dissipation. In section 2 and 3 we 
briefly describe the equations and the numerical method used and we present the simulations that were 
performed. In section 4 we present an overall picture of the dissipation fields through the time evolution 
of their average values, their pdfs and their spectra. We also provide a qualitative view of the dissipation 
field in physical space through color maps. Section 5 deals with the extreme dissipative events, it begins 
with a discussion of the structure functions of the velocity and magnetic fields and concludes with results 
from the statistical analysis of the geometrical and dynamical properties of structures of high dissipation. 
We discuss and conclude our results in section 6. 


2 The equations 


2.1 Ambipolar drift 

In a partly ionized fluid, the time-dependent evolution of both the neutral and the ionized fluids should 


in principle be followed. However, in circumstances that we will make explicit below (see subsection 2.4), 
we can adopt the strong coupling approximation. In this approximation, we neglect the inertia, pressure 
and viscosity of the ions in the ion momentum evolution and we are left with the balance between the 
ion-neutral drag and the Lorentz force 


0 / PzPn(Ui U n ) — J X B 


( 1 ) 


where the current density 

J = i-(VxE), 

Pi and p n are the ion and neutral mass density (pi <C p n ), U* and U n are the ions and neutrals respective 
velocities, and where 7 is the coefficient of ion-neutral drag 


(crv)in 
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with rrii and p the ions and neutrals mass per particle and {av)i n the ion-neutrals collision rate. Assuming 
that 77-He = 0.2nH 2? we find p = 2.33 m v for molecular gas, where m p is the mass of t he proton . In diffuse 


1983 


clouds, the average mass per ion is m* = 12m p as the dominant ion is C + . Following Draine et al. 
we take ( crv)i n = 1.9 x 10 - 9 cm 3 s -1 and we arrive at 7 = 6.7 x 10 13 cm 3 s -1 g -1 . 

Within these approximations, the above balance ([Tj) expresses the ion-neutral drift velocity as a 
function of the magnetic field. When this is plugged into the induction equation 

a t B = Vx (Ui x B) +ryV 2 B 


one recovers the ambipolar diffusion term: 


V x 


7 PiPn 


-(J x B) x B 


[cf. [Balbus an d Terquem 2001], which can be developed into 

V x 




yPiPn 7 PiPn 


from which properly speaking onl y the first te rm takes the form of a diffusion, with diffusion coefficient 


Aad — 


Brandenburg and Zweibel 


1994 


but the qualificative is usually retained for the whole 


term. In particular, Brandenburg and Zweibel [1994] recognized that the second term steepens the 


magnetic field profile near magnetic nulls. It should finally be noted that AD itself is not per se able to 
reconnect the field fines: this requires Ohmic diffusion. 


2.2 Incompressible MHD 

We now take xxo a typical velocity and lo a typical length scale as unit velocity and unit length to 
normalize our equations. We also define to = lo/uo as the unit of time. We write the non-dimensional 
velocity u = U c dm/^o where U c dm is the center of mass velocity 


U c dm — 


PiUj + Pn^n 


with p rsj p n the total density of the gas, p — pi + p n - We write the non-dimensional Alfven velocity 
b = B/a/47tp/xxo. The non-dimensional current is simply j = V X b where V is now understood as 
derivatives in coordinates in units of lo: V ZoV. 

In the diffuse ISM, the sonic Mach number M s — uo/c s (where c s is the s peed of sound) as well 
as the Alfven Mach number M a — l/|b| take values in the range 10 -1 — 10 Elmegreen and Scalo 


2004 . This wide range of values suggests that although most of ISM turbulence is highly compressible 
incompressible turbulence is not irrelevant since in a weakly compressible flow the density fluctuations 
A p/p ~ Mg are of the order of the square of the sonic Mach number. Hence a turbulent flow with 
Ms < 0.3 can be adequately described by the incompressible equations. For example, the studies of 
|Branden burg and Zweibel [1995], 


Zweib el] [2002] or | Godard et al.| [2009] on turbulent dissipation with 


AD were all based on the incompressible equations of motion. 

We use the equations in Balbus and Terquem [2001 and the above notations to derive the equations 
of incompressible, viscous, resistive, AD MHD : 


d t u + (u • V) u = — V p + j x b + Re 1 V 2 u 

dtb = V x (u x b) + Re~ x V x ( (j x b) x b) + Re^ 1 V 2 b 


where u and b satisfy V • u = 0 and V • b = 0 and p m. P/(uqp) is the non-dimensional pressure, with 
P the actual thermal pressure. Equations § are parameterized by three non-dimensional numbers Re , 
Re a and Rem , for which we now give estimates. 
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2.3 Reynolds numbers 

Firstly, the kinetic Reynolds number Re = uolo/v, where v is the molecular viscosity of the fluid, expresses 
the relative importance of inertial terms in comparison to the viscous term. In the neutral ISM, assuming 
that the most significant contribution to viscosity is given by H 2 collisions, we have v ~ |Ah 2 c s where 
Ah 2 is the mean free path of the H 2 molecule and c s = (VksT / p) 1 ^ 2 is the isentropic sound speed, with 
r ~ 5/3 the ratio of specific heats, and ks the Boltzmann constant. The mean free path is given by 
Ah 2 ~ (mi 2 cr h 2 ) -1 where tih 2 is the number density of H 2 and cth 2 = 3 x 10“ 15 cm 2 is an estimate of 
the cross section of H 2 collisions [Monchick an d Schaefer, 1980 . For molecular gas nu 2 — 0.5riH where 
riH is the hydrogen nuclei density. Under these assumptions, the kinetic Reynolds number is of the order 


Re 


18 x 10 ' (srS^) (tiSA) 


Iq 


10 pc ) \ 100 K 


Elmegreen and Scalo 2004 quote typical values of the kinetic Reynolds number in the cold ISM ranging 


from 10 ° to 10 Y . 

Secondly, Re m = uolo/rj is the magnetic Reynolds number, where 77 is the resistivity. The magnetic 
Reynolds number expresses the relative importance of advection in comparison to Ohmic diffusion in the 
dynamics of the magnetic field. In a system with a large value of the magnetic Reynolds number, the 
dynamics of the magnetic field is dominated by advection and stretching. The value of the resistivity is 
given by 

■n = 234 (—^ T 1/2 c 


[Balb us and Terquem 2001] , where n is the total number density and n e is the electron density. The 
order of magnitude of the magnetic Reynolds number is 


Rem = 2.2 x 10 1 


f l ° ^ ( u ° 

\ 10 pc/ v 1 km s _1 / 
T 


io- 




100 K 


km s _1 
- 1/2 


where we assumed n = 0 . 6 nH for molecular gas. 

Lastly, the AD Reynolds number Re a helps to measure the ratio of the ambipolar to advective 
electromotive forces in the induction equation: 


R&a — , 5 

ta 


ta 


If Pi 


(3) 


where t a can be recognized in equation |l]) as the ion-neutral friction time scale. The quantities Re a 
and l a should not be confused with their more usual definitions Rad and £ad as introduced by |ZweibeI| 
and Brandenburg [1997], for example. For instance, the usual values depend on the r.m.s. velocity and 
magnetic field, whereas our definition encompasses only the ion-drift time: see the subsection |2.4| for 
more details. We can also write it as a ratio of length scales 

v _lo 
R&a — , 


where we define 


la — t a U 0 


(4) 


which gives a typical length scale for ion-neutral decoupling. The AD Reynolds number is an increasing 
function of the ionization fraction x — pi/p. Using C + as the dominant ion of the ISM, we find 


Re a = 4.9 x 10 3 



( n c+ \ 

/ n H W \~ 1 

g 

1 

o 

V100 cm -3 / VI km s -1 / 


In the ISM, we have Re a <f?e< Re m which suggests the ordering l a l v l v for the ambipolar, 
viscous and resistive dissipation scales. However, our finite computing power does not allow much 
dynamics of scales and here we can only afford l a > Id where Id = l u ~ 1-q is a single dissipative scale. 
In this study the magnetic Prandtl number Pr m = v/r\ is therefore taken equal to unity, so that the 
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hydrodynamic and magnetic Reynolds numbers are equal. This choice results in a single dissipative range 
of scales for both the velocity and the magnetic fields, a fact which simplifies the analysis considerably. 

Although this numerical study is confined to Reynolds numbers that are several orders of magnitude 
lower than those found in the interstellar medium, it is relevant to ISM physics in the sense that it allows 
a detailed quantitative study of the dissipation field as well as the relative importance of the different 
types of turbulent dissipation, which can all be important as a heating source for ISM chemistry. 


2.4 Lengths scales associated with AD 

We outline here various scales introduced by AD physics. The case of C-shocks allows to clearly separate 
them. We have already introduced the length scale l a = uot a which corresponds to the re-coupling length 
scale between ions and neutrals. It is the scale of variation of the neutral’s velocity in a C-shock with 
entrance velocity uq , which therefore has a length on the order of l a [Flower and Pineau des Forets 


1995 


Under the strong coupling approximation, the scale of variation of the ions velocity in the same C- 
shock is l a i = uota/Ma where A4 a is the transverse Alfvenic Mach number of the shock. For shocks with 
a large Alfvenic Mach number, this length-scale is significantly smaller than l a and the structure of the 
shock consists of a front in the ions velocity followed by a smoother transition for the neutral velocity 
[see Fig. 1 of Li et al. 2006, for example]. 

In the case of typical ISM turbulence, though, M a is of order one, and both length scales do not differ 
significantly. In fact, jZweibel and Brandenburg 11997 constructed the AD diffusion Reynolds number of 
eddies of length scale i and velocity U bas ed on the AD diffusion coefficient : Rad(£) = IU / Aad = ^/4d 
where l ad — Ut a /Ma with M a — ET/ca |j Zweibel and Brandenburg 1997 then argue that only eddies 


of length scales below 4 d should be affected by AD. We prefer to get a similar estimate by comparing the 
Fourier amplitudes of the AD e.m.f. (Re” 1 (j x b) x b 4 Re~ 1 kb 2> ) and the inertial e.m.f. (uxb4 ub) 


the the induction equation Wave numbers above the critical wave number 

k a = Re a \J(y 2 )/^) 2 ) 


(5) 


should be AD dominated. We hence define £ a = 27t /k a accordingly, as the length scale below which AD 
should be effective. Note that 4 d and lo£ a differ only by a factor of 27r. 

Similarly, we can estimate the length scale below which the strong coupling approximation breaks 
down by comparing the magnitude of the neglected inertial term DtpiXJi to the coupling term p n (U n — 
U i)/t a . Assuming that U n , U* and U n — U* all share the typical magnitude uq, we then get a critical 
wavenumber /ctwo-fluids — pn/Pi/la above which the strong coupling approximation fails and the two- 
fluids approximation is needed. Provided pi/p n is small (it is typically lower than 10 -3 if the main 
charges are C + ions), the strong coupling approximation breaks down at scales much smaller than the 
typical AD diffusion scale. Other authors [Oishi and Mac Low) 2006, Padoan et al. 2000 have claimed 
that the strong coupling approximation breaks down as soon as £ < £ad or Rad < 1, where the ions 
inertia does not appear explicitly. But Fig. 1 of Li et al. 2006 shows a C-shock computed with the 
two-fluid approximation (solid) compared to an analytical solution (dashed) using the strong coupling 
approximation, and the agreement is perfect. We hence believe that the strong coupling approximation 
is a very good one in the low ionized ISM where pi/p n <C 1, even in cases where Rad > 1 . In particular, 
^two-fluids is at least a few ten times above the largest wave number in all our AD simulations, which 
amply justifies our use of the strong coupling approximation. 

Finally, the observed emission of the ISM depends on the chemical and thermal state of the gas, 
which are strongly linked to the heating. The scale at which the heating takes place may not necessarily 
be directly connected to the scale £ a where AD undergoes a change of dynamical regime. In fact we will 
see that it is not the case in the present paper, and we will be forced to introduce yet an other length 
scale £* for the typical thickness of sheets of strong AD heating. 


x For example, with i m l a , we find that the AD Reynolds number of a C-shock is either Rad = 1 or Rad = At^ depending 
on whether we compute it with the post-shock or the pre-shock magnetic field strength. 
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3 The simulations 

3.1 Method 


Many compressible methods have been devised to treat AD in the strong coupling or in the two-fluids 
approximations: see Masson et al. 2012] and references therein. Here, we solve the strong-coupling 


incompressible equations ([2| in 3D using a spectral method for various values of the parameters and 
various initial conditions. Our spectral code ANb0is fully de-aliased by use of the phase-shift method of 
Patterson and Orszag 1971 and uses polyhedral truncation [Canuto et al. 1988| . Polyhedral truncation 
considers only these wave-vectors for which the sum of any two of their components does not exceed 
2N/3. Similarly with the widely used two-thirds rule, but contrary to an isotropic spherical truncation, 
this truncation scheme does not possess a sharp limit in wavenumber space. Polyhedral truncation allows 
us to keep 55 percent of the modes active, in comparison to 33 percent for the standard two-thirds rule, 
resulting in a more accurate description of the small scales. The Fourier transforms are computed with 
FFTW with single precision accuracy, and a standard fourth-order Runge-Kutta method is used for time 
integration. We checked that our code gives the correct solution for Alfven waves damped by AD. As a 
resolution check for all simulations, we also looked for bumps in the kinetic and magnetic energy spectra 
near the truncation limit. In the case of the kinetic energy spectra we found bumps no larger than 15 %, 
whereas no bumps were present in the magnetic energy spectra. Note that we do not include a driving 
force in equations our simulations are freely decaying. The MHD simulations with 512 3 resolution 
take about 5000 CPU hours until the peak of dissipation but the equivalent AD-MHD simulations require 
ten times more CPU time due to the more stringent time-step requirement. 

In the spectral simulations performed, the velocity and magnetic fields are defined on a regular 
Cartesian grid of points, while boundary conditions are periodic in all directions. Note the total length 
of the computational domain is 2tt and the smallest non-zero wave-vector has a norm of one. 


3.2 Initial conditions 


We use two types of initial conditions, corresponding to two different situations for the magnetic and 
cross-helicities. 

In the first case, the three lowest non-zero wave numbers of both the velocity and the magnetic field 
are initially loaded with a superposition of different Arnol’d-Beltrami-Childress [ABC, see Dom bre et al.[ 
1 986] flows 


(u x ,u y ,u z ) = (Asm(kz) + Ccos(ky), B sm(kx) 
+ A cos (kz ), C sin (ky) + B cos (kx)). 


(6) 


Different values of the coefficients A, F>, C are chosen for the first three wave numbers from a uniform 
random number generator. In higher wave numbers a random field with energy spectrum 

E{k) = C E k~ 3 exp (-2 (fe/fe c ) 2 ) , fe c = 3 (7) 

is superposed. The phases are chosen from a uniform random number generator with the same seed for 
all simulations. 

In the second case, the large scale initial flow is the Orszag-Tang (OT) vortex 


(■ u x ,u y ,u z ) = (—2siny, 2sinx, 0) 

(b X: b y ,b z ) = (~2sm(2y) + sinz, (8) 

2 sin(x) + sin z, sin x + sin y) 

and in higher wave numbers a random velocity field with the same properties as above is superposed. In 
order to keep the initial value of magnetic helicity close to zero, no random magnetic field is added to 
the OT initial condition, in contrast to the ABC initial condition. 

The compressive components of the initial velocity and magnetic fields are subtracted so that the 
initial condition is purely solenoidal. In all cases, the constant Ce in equation 0 is chosen such that 

2 http://www.Ira.ens.fr/~giorgos/ank 
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# 

N 

L 

A 

i d 

Re x 

S 

03 

0$ 

II 

03 

0$ 


Initial condition 

1 

128 

2.60 

1.30 

0.0607 

189 

219 

- 

ABC 

2 

128 

2.59 

1.31 

0.0605 

210 

219 

- 

OT 

3 

128 

2.70 

1.44 

0.0390 

217 

219 

100 

ABC 

4 

128 

2.72 

1.44 

0.0381 

228 

219 

100 

OT 

5 

256 

2.46 

0.94 

0.0375 

353 

551 

- 

ABC 

6 

256 

2.46 

0.91 

0.0389 

377 

551 

- 

OT 

7 

256 

2.62 

1.08 

0.0185 

412 

551 

100 

ABC 

8 

256 

2.66 

1.09 

0.0196 

444 

551 

100 

OT 

9 

512 

2.28 

0.67 

0.0237 

591 

1374 

- 

ABC 

10 

512 

2.12 

0.58 

0.0245 

604 

1374 

- 

OT 

11 

512 

2.49 

0.83 

0.0091 

756 

1374 

100 

ABC 

12 

512 

2.36 

0.75 

0.0095 

750 

1374 

100 

OT 

13 

512 

1.84 

0.58 

0.0074 

640 

1374 

10 

ABC 

14 

512 

2.80 

1.00 

0.0084 

927 

1374 

10 

OT 


Table 1: Parameters of the simulations. N: linear resolution, L : integral length scale at the peak of 
dissipation (pd) , A: Taylor microscale (pd), dissipative scale (pd), Re\: Taylor microscale Reynolds 
number UXRe (pd), Re\ kinetic Reynolds number, Re m : magnetic Reynolds number, Re a : AD Reynolds 
number. 


(u 2 ) = <b 2 > = 1, so that we start from equipartition between kinetic and magnetic energy. The energy 
of the initial condition fields is concentrated on large scales k < k c due to the exponential cutoff in 0- 
In the case of the ABC initial condition, the non-dimensional cross-helicity 

H = 2(u • b) 

V< u2 X b2 > 

is ~ 2 x 10 —3 , corresponding to a low initial correlation between the velocity field and the magnetic field. 
The mean magnetic helicity 

Hm = (a • b) 

where a is the vector potential with b = V x a, is considerable, ~ 0.2. In the case of the OT initial 
condition the non-dimensional cross-helicity is ~ 0.1 while the mean magnetic helicity is almost zero, 
~ 1 x 10 -9 . Thus these two different initial conditions represent evolution under different constraints: 
in the ABC case, low cross-helicity and sizable magnetic helicity whereas in the OT case sizable cross- 
helicity but low magnetic helicity. This fact is important because in the ideal MHD limit (inviscid and 
non-resistive) the energy, cross-helicity and magnetic helicity are all conserved during the evolution. If 
AD is included in the ideal MHD equations, the conservation of magnetic helicity remains while energy 
and cross-helicity conservation are broken. This is a consequence of the form of the AD term in the 
induction equation, which takes the form of an advection term 

V x (ud x b) , with u d = Re~ x { j x b) 

the non-dimensional ion-neutral drift velocity. This form also implies that although AD is a dissipative 
process, it conserves magnetic flux and is thus unable to reconnect field lines. 


3.3 Parameters 


The parameters of the simulations performed are shown in table [I] Throughout this paper, we focus 
mainly on the analysis of the OT initial condition with Re a = 10,100, oo, and we discuss the differences 
with respect to the ABC initial conditions only when they arise. The Taylor microscale Reynolds number 
Re\ is defined as Re\ = UXRe where U = \J (u 2 ) is the r.m.s. velocity and 


A = 2tv 


( r e ( fc ) 

v/q 00 k 2 e(k) dk 
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is the Taylor microscale, with J 0 °° e(k) dk = |(u 2 + b 2 ) the total energy and e(k) the total energy 
spectrum. The value of Re\ is given at the peak of viscous plus Ohmic dissipation where 


(e> = ( £ «> + fa) 

_ T~) — 1 • 2 

So R&m J 


£v = Y! ( diU i + d J u i ) 2 

i,j —1 


The time when the peak of dissipation occurs is appropriate for analysis because for a given integral 
length scale 


r _ o Jo°° k le ( fe ) dk 
f 0 °°e(k)dk 


(9) 


and dissipative scale (assuming Kolmogorov scaling) 


ld = 



( 10 ) 


the scale separation L/ld between the energy-containing scales and the dissipative scales is maximum. 
Another desirable property at the peak of dissipation is quasi-stationarity, due to the time derivative of 
the dissipation rate which cancels at the peak, by definition. In all the MHD simulations (with Re— 0), 
two snapshots of the fields were recorded for analysis: one at the peak of dissipation and a second one 
one eddy turnover time later. We define the macroscopic eddy turnover time as 

T= (11) 

where the one-directional r.m.s. velocity U/V?> and the integral length scale L were both computed at 
dissipation peak to estimate when to output the next snapshot. For the AD simulations (with Re >0), 
we could not afford to compute beyond the dissipation peak. 


3.4 Power-spectra 

The kinetic and magnetic energy spectra of the high resolution OT runs 10-12-14 are shown in Figure 
[l] at the temporal peak of dissipation. The spectra are compensated by the Kolmogorov law k ~ 5 / 3 and 
normalized by U 2 . The extent of the inertial range, as defined by the portion of the spectra that has 
slope —5/3 is very limited, especially in the cases with AD. In the same Figure we show the limits of the 
inertial and dissipation ranges assumed for the analysis of section [572] they are taken from |Uritsky et ah] 
[2010 as [0.21,1.3] and [0.025,0.18] respectively (in units of Zo). 

The kinetic and the magnetic energy appear to remain in approximate equipartition across all scales 
except for the smallest scales in the AD runs where magnetic energy dominates the kinetic energy (the 
tick-marks on the vertical lines can guide the eye to estimate the relative position of the curves between 
the upper and the bottom panel). 

The vertical dash-dotted lines with square symbols correspond to the wave-number k a defined in 
equation for the cases of runs 12 and 14. Surprisingly, departures from MHD spectra start at about 
the same wave number for all AD runs (in the range k ~ 5 — 8 for both the kinetic energy spectrum and 
the magnetic one): this hints at the fact that £ a — 2n/k a is not the proper scale to assess the dynamical 
importance of AD in our simulations. In particular, dynamics can be affected at scales much larger than 
that in run 12. Although the difference between the pure MHD and AD MHD spectra is modest, there 
is a clear tendency for AD to flatten the energy spectra, especially for the magnetic energy. This is in 
line with the idea by Brandenburg and Zweibel 1994 that AD diffuses magnetic fields on the one hand, 
but on the other hand helps to build sharper magnetic structures in specific places. 
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Figure 1: Compensated kinetic (top) and magnetic (bottom) energy spectra for OT runs 10, 12 and 14 at 
the temporal peak of dissipation. The assumed limits of the inertial and dissipation ranges are shown in 
dashed vertical lines. The vertical blue and red dashed-dotted lines with square symbols correspond to the 
expected AD critical wavenumber k a (see equation <§) in runs 12 and 14, respectively. 
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Figure 2: Time evolution of volume integrated dissipation rates for the OT runs 10, 12 and 14. Solid lines 
show the Ohmic plus viscous dissipation which we use to define peak dissipation. 


4 The dissipation field 


4.1 Total dissipation 

We present in Figure [2]the time evolution of volume integrated dissipation rates. The viscous and Ohmic 
dissipation rates follow each other closely and we don’t separate their respective contributions in this 
figure. Most of the time, the total dissipation rates due to viscosity, resistivity and AD are of comparable 
magnitude. But the total AD dissipation rate is seen to peak before the Ohmic plus viscous dissipation 
rate, especially at low values of Re a . This makes our choice of the temporal peak of Ohmic plus viscous 
dissipation more appropriate to avoid the initial transient spike of AD dissipation (this spike is even more 
pronounced in the ABC case run 13). 

We note that Ohmic dissipation is not enhanced by the presence of AD. On the contrary, the peak 
value of the Ohmic plus viscous dissipation decreases as Re a is decreased. Even though [Brandenburg and| 
ZweibeT jl994]’s idea that AD sharpens magnetic structures at small scales is valid in our simulations, AD 
also smooths the fields at intermediate scales and the net effect on the global Ohmic heating is to decrease 
it. However, this may be due to the finite dynamical range in our simulations. Higher Rem Reynolds 
number simulations, if they yield enhanced magnetic power in a more extended range at small scales, 
could result in a globally enhanced rate of reconnection, in agreement with |Zweibel and Brandenburg| 
1997 ’s model. 
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Figure 3: Log-normal core and power-law tail fit for the pdf of the total dissipation (Run 12 AD - OT). 
The core follows a log-normal distribution with mean Hi ~ —4.27 and standard deviation cp ~ 1.03 while the 
tail follows a power-law with exponent —r ~ —2.61. Vertical lines show the mean value (red) and thresholds 
located at 1 (black), 2 (green) and 3 (blue) standard deviations above the mean value. 

4.2 Probability distribution function 

The probability density function (pdf) of the total dissipation rate 


St — s 0 T s v T £ a 


where 

£a = R( J a, 1 (j X b) 2 

is shown in Figure [3] for the high-resolution run 12. The core of the pdf is very close to the log-normal 
distribution 

rr-y / \ ( (lll£ f -/ii) 2 \ 

V c {£t) oc exp - ^2 - — J 

with mean m ss —4.27 and standard deviation cp ~ 1.03, while the tail of the distribution can be fitted 
by a power-law 

Vt(st) oc £t T 

with exponent r ~ 2.61. This power-law is one of the signatures of intermittency of dissipation |Frisch[ 
1995]. For still higher values of the total dissipation the pdf has an exponential cut-off, although these 
high dissipation values are close to the sampling limit. In the analysis of the next section (extraction of 
structures of high dissipation) effectively only the power-law range of the distribution is sampled. 

In Figure [4] we present the cumulative probability density function of the total dissipation for run 12. 
It flattens out before reaching unity, which shows that high values of the dissipation are concentrated in 
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Figure 4: Cumulative probability density function of the total dissipation for run 12 (AD - OT). 


a small volume subset of the spatial domain. The analysis of the next section concerns events that take 
place in this high-dissipation plateau. 


4.3 Power-spectrum 

We now investigate the distribution of the energy dissipation with respect to spatial scale. This can be 
discerned thanks to the power spectra of the dissipation fields, displayed in figure [5] The AD heating 
peaks at larger scales in comparison to Ohmic and viscous dissipation. The scale of this peak £* is 
actually only a few times smaller than the integral length scale. At this range of scales AD dissipation 
is much more important than Ohmic and viscous dissipation. This suggests that the heating due to AD 
has a characteristic length scale i* a , which can be much larger than the dimensional estimate i a — Zn/ka 
§ for run 11. This length scale is relevant to the heating, and hence may manifest itself in structures 
revealed by chemical tracers of the warm chemistry of the ISM. 

The AD dissipation term, i?ey 1 ( j X b) 2 , is proportional to the square of the Lorentz force. It is 
hence interesting to look at the influence of AD on the characteristics of the Lorentz force, in comparison 
with the pure MHD case. Figure [6] shows the spectrum of the Lorentz force for OT runs 10, 12 and 
14. The inclusion of AD has a significant effect on the total power of the Lorentz force, reducing it 
importantly especially in the dissipative range. By contrast, as seen on figure^ AD results in a deficit 
in the magnetic energy spectrum (and hence the spectrum of the current vector) which is much smaller 
in comparison to the deficit of the cross-product of these two vectors (the Lorentz force), and is only 
present on intermediate scales. 

This is explained if AD has the effect of aligning the vectors j and b [as was also found in the 
simulations of Brandenburg et ah] 1995 , with a stronger effect at small scales. To put it in an other 


way, AD leads the magnetic field at small scales closer to a Lorentz force free configuration, where the 
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Figure 5: Compensated dissipation spectra for run 12 (AD-OT, Re a = 100). Blue solid line: AD dissipation, 
red dashed line: Ohmic dissipation, green dotted line: viscous dissipation. We plot ke e (k ) in a log-lin plot, 
so that the area under the curve over any interval shows directly the amount of power inside this interval. 
We mark the position of the maximum value of ke eal at k = k* = 2 tt/£*, and the position of k = k a . 
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Figure 6: Power spectra of j x b for high resolution runs 10,12 and 14. The field becomes force-free at small 
scales when the strength of the AD is increased. 


feedback of the magnetic field evolution on the velocity field dynamics is weaker than for MHD. This 
was also found in the simulations by Brandenburg and Zweibel [1995] . 

Here we attempt to trace this tendency back to the evolution equations. We write jj_ the component 
of the current vector j: the double cross-product (j x b) x b can then be simply written b 2 jj_. This allows 
to write 

(ftj)withAD = (^j) without AD + V X [V X jx)] (12) 

which shows that in regions where b 2 is smooth enough, the effect of AD is to diffuse out the component 
of the current perpendicular to b, and it does so faster at small scales like any diffusion process. Hence 
AD brings the field closer to a force-free state, more efficiently at small scales, except perhaps at the 
smallest scales where b 2 varies and the behavior of equation (12) is less easy to predict. 


4.4 Spatial structure 

Next, we consider qualitatively the different contributions to the bulk of the dissipation field. For this 
purpose, each different mechanism of dissipation (Ohmic, viscous and ambipolar diffusion) is assigned to 
a color channel: Ohmic dissipation is assigned to the red channel, viscous dissipation to the green channel 
and AD dissipation to the blue channel. To emphasize the structures in the bulk of the dissipation we 
first compute the total dissipation value £i below which 10% of the heating occurs and the value e u below 
which 90% of the dissipation occurs. We discard the pixels with total dissipation et < £i, we saturate 
the intensity of the pixels with e t > £ u (while keeping their intrinsic color) and we apply a logarithmic 
scaling for the intensity in between these two thresholds. The color of each channel is hence given by the 
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Figure 7: Color maps of a slice of the dissipation fields, run 10 (MHD - OT). Red: Ohmic dissipation, green: 
viscous dissipation. All snapshots are taken at the peak of dissipation. 


ratio of each type of dissipation to the total dissipation 


with the intensity I given by 


I = 


Red = — I 

St 

Green = — I 
S t 

Blue = — I 

s t 


0 

l°gQt) — log(ez) 

log(e u )-log(ez) 

1 


if s t < si 
if si <s t < s u 
if e u < s t 


The color maps of a slice through the dissipation fields are shown in figures [7|[9| for all high-resolution 
runs with the ABC initial conditions. 

In the pure MHD case (Figure [7]), viscous and Ohmic dissipation are in general concentrated on 
thin sheets: a slice by slice inspection of the full cube reveals continuously evolving filaments at the 
intersection between the sheets and the plane of the slice. The length of the sheets is comparable to the 
integral length scale 0 while their thickness is comparable to the dissipation scale |To| . 

The case of AD MHD (figures [8 1 9 [ ) is similar in the sense that viscous and Ohmic dissipation are again 
concentrated on thin current sheets, although these sheets are fewer in number (more connected), and 
the voids of low dissipation between them are smaller. AD dissipation is significantly more diffuse than 
both Ohmic and viscous dissipation, and is concentrated on much thicker structures. The thickness of 
the AD dissipation structures seems to coincide with the AD heating length £* measured on the power 
spectra. In some cases the structures of strong AD dissipation are seen to surround the sheets of Ohmic 
or viscous dissipation: AD sandwiches Ohmic dissipation, much like in the Brandenburg and Zwe ibel 
[1994] picture. This is probably also seen on the left panel of Figure 7 of Padoan et al.| [2000] where the 
structures of AD dissipation often go in pairs, except Ohmic dissipation is absent from their simulations, 
so reconnection proceeds only through numerical truncation errors. Between figure [7] and Figure [§] there 
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Figure 8: Same as Figure [7| Run 12, (AD - OT, Re a = 100) with blue: ambipolar diffusion heating. 



Figure 9: Same as Figure [7[ Run 14, (AD - OT, Re a = 10) 
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is little difference in the spatial structure of the dissipation field, with similar size, shape and position 
for each sheets of dissipation. 

The above color maps reveal only a small fraction of the pixels have mixed colors (such as cyan, 
yellow or purple). This suggests that the dissipative structures of different nature (Ohmic, viscous or 
AD) are well separated. 


5 Intermittency and structures of high dissipation 

5.1 Structure functions 

Figures [7]|9] suggest that the dissipation field is not smoothly distributed in space, but is characterized 
by a high degree of intermittency with regions of extreme dissipation values alternating with relatively 
quiescent regions. The intermittent distribution of the energy dissipation rate is expected to have an 
impact on the structure functions of the velocity and the magnetic field. A longitudinal velocity structure 
function of order p is the p -th moment of the longitudinal increment of the velocity field 


Sp(r) = ((<5«||(r)) I> ), <5w||(r) = (u(x + r) - u(x)) • r 

where r is the unit vector in the direction of r. The definition for t he st ructure function of the magnetic 
held Sp(r) is completely analogous. According to standard Kolmogorov 1941 phenomenology (hereafter 


K41), in the inertial range the structure functions exhibit power-law scaling 


Sp (r) oc r^ p , 


Id « r < L 


where the exponents (p vary linearly with p, = p/3. Intermittency considerations [Frisch 1995] lead 
to deviations from the linear K41 prediction. A particularly successful model of intermittency that was 
introduced for hydrodynamic turbulence by She and Leveque 1994 and generalized for MHD turbulence 
by Politano and Pouquet [1995] predicts 


£° L (9,C) =-(!-- 




9 \ 9y 

where g is the inverse of the inertial range scaling exponent of the velocity increment 

5u\\ (r) oc r 1 ^ 9 

and C is the co-dimension of the dissipative structures, C — 2 for filaments and C — 1 for sheets in 
three space dimensions. K41 phenomenology predicts g — 3, while Iroshnikov-Kraichnan (IK) MHD 
phenomenology |Iroshn iko~v| 1964 Kraichnan 1965] predicts g — 4. 

In order to estimate the structure function exponents from numerical data, we used the extended self¬ 
similarity (ESS) method introduced by Benzi et al. 1993 . The exponents were estimated by computing 
the logarithmic slope 

= d log .Sp (r) 
p d log S'* (r) 

where S*(r) is a structure function whose scaling behavior is known from theory. In the case of MHD 
turbulence, S*(r) is given by the two functions 

St(r) = (5^(<5z ± ) 2 } 

where Sz± is the increment of the Elsasser fields z± = u ± b. Starting from the MHD equations and 
assuming statistical homogeneity, isotropy and stationarity, one can derive analytically 


- -^> 


Id < r < L 


(13) 


where e ± are the dissipation rates of (z^) 2 Politano and Pouquet] 1998]. The linear scaling of Sf(r) 


the inertial range is confirmed approximately by the flattening seen in Figure [lO] which shows compensated 
plots. Although the derivation of the law (13) is not proven in the case of AD MHD, the linear scaling 
of S±(r) with r is not further from linearity in comparison to pure MHD. 
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Figure 10: Compensated plot S±(r)/r for the OT runs 10,12 and 14. 
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The velocity and magnetic field structure function exponents up to order 8, calculated for high 
resolution runs 9-14 using extended self-similarity are shown in figures pTp4| We see a departure from 
the linear prediction of K41, a sign of intermittency, especially in the case of the magnetic field. The 
structure function exponents follow closely, but not exactly, the predictions of the generalized She &; 
Leveque model with a scaling parameter g — 3 and a co-dimension for the structures of high dissipation 
around C — 1 (with exceptions at C = 2 and some below (7=1). This fact suggests that in the inertial 
range Su(r) and Sb(r) are proportional to r 1//3 , which is the prediction of K41 phenomenology, rather than 
r 1 / 4 as predicted by IK phenomenology. The value C = 1 of the co-dimension suggests that the structures 
of high dissipation in the inertial range are sheet-like, in accordance with figures |7|9| and the analysis 
of section 5.2 The velocity exponents for runs 9, 10 and 12 hint towards (7 = 2 (filaments) and the 


magnetic exponents for the OT runs exhibit a greater degree of intermittency. This indicates that even 
though the exponents appear to follow generalized She & Leveque models, the model’s phenomenology 
probably does not subtend the dissipation in our simulations. 

The effect of AD on intermittency is not easy to discern from these scaling exponents. In the case 
of the ABC initial condition, the deviation from the K41 values is larger for both the velocity and the 
magnetic field. The situation is reversed in the case of the OT initial condition, where the pure MHD 
fields appear to be more intermittent. 


5.2 Extraction of structures 


Following Uritsky et al. polo] (hereafter UR10), we implemented an algorithm for the extraction of 
structures of high dissipation. In the problem considered, there are three types of local dissipation rates: 
the viscous dissipation rate, 


the Ohmic dissipation rate 
and the AD dissipation rate 


£ v = ( diU i + d i Ui ) 2 

i,j=l 

£o = -Re” 1 j 2 = fie” 1 (V x b) 2 


£ a = Re a 1 (j x b) 2 . 

A structure of high dissipation is defined as a connected set of points x where 

e(x) > (e) + ja e , j =1,2 or 3. 


(14) 

(15) 

(16) 

(17) 


In the above relation, e can be any of the three dissipation rates or the total dissipation rate St — 
s v + So + s a (or e t — e v + So in the case AD is absent). (e) is the spatial mean value of the dissipation 
rate and a £ its standard deviation. For example, the three thresholds we use on total dissipation for run 
12 are displayed on figures [3] and [4] over the PDF and the CDF of the total dissipation. 

The algorithm is capable of isolating the structures of high dissipation so that a statistical analysis of 
their geometric and dynamical characteristics can be performed. The extracted structures are generally 
sheet-like, as can be seen in figure |15[ where all structures extracted from run 12 (AD - OT) whose 
characteristic linear sizes Li (see below) lie in the inertial range are shown. An example of a more 
complex structure can be seen in Figure |16| This structure is one of the largest extracted from this 
dataset. It is sheet-like, with its length larger than the integral length scale § and its thickness is 
comparable to the dissipation scale (10). It is overall characterized by a high degree of geometrical 
complexity. 
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Figure 11: ESS velocity field structure function exponents for ABC runs 9-11-13. 
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Figure 12: ESS velocity field structure function exponents for OT runs 10-12-14. 
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Figure 13: ESS magnetic field structure functions for ABC runs 9-11-13. 


22 







— 

K41 

- 

- 

(p SL (3, 2) 

,,, 

■ ■ 

Cp SL (3, 1) 

☆ 


C p GSL (4,1) 


it 

( GSL (4, 2) 

X 

X 

Run 10 (MHD 

☆ 

☆ 

Run 12 (AD O' 

• 

• 

Run 14 (AD O' 


Figure 14: ESS magnetic field structure functions exponents for OT runs 10-12-14. 
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Figure 15: Inertial range structures extracted from the dataset corresponding to the peak of dissipation of 
Run 12 (AD - OT). They are defined as connected sets of points having values of the total dissipation three 
standard deviations above the mean value. Each little sphere has a 2-pixels diameter, ie: about the size of 
the viscous (or equivalently Ohmic) dissipation length. 



Figure 16: One of the largest structures extracted from run 12 
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After the extraction of the structures, the following quantities were computed for statistical analysis 


Li = 5 max |r m - r/| (18a) 

m,ZGA^ 

L{ x , v ,z}i = 5 max \{x,y,z} m ~ {x,y,z}i\ (18b) 

m,tG 

Ri = \[PXx~XL zi (18c) 

Vi = S 3 1 ( 18d ) 

m£Ai 

Ai = 5 2 1 ( 18e ) 

mGAj 
M (m) (£. A-j 

Pi = 5 3 J2 £ ( r m) (18f) 

mGAj 


In the above definitions, A* is the i-th structure, m its ra-th point, M(k) the set of 26 neighbors of the 
point m and 5 = 2ir/N is the grid spacing. Li is the characteristic linear scale of the structure, L{ XjVjZ yi 
is the linear scale of its projection on the three axes, Ri is the characteristic linear scale of the smallest 
volume embedding the whole structure, Vi is its volume, Ai its surface and Pi the volume-integrated 
dissipation rate. In the next section, all statistical quantities are computed on the sample of different 
extracted structures. Note that the definition of Ri makes it dependent of the orientation of the structure 
(its value can change by tilting slightly the axis of the domain). 

In order to validate our extraction procedure, we implemented it in two different ways. First, we 
implemented the same algorithm as UR10 (ie: recursive, using breadth-first search as explained in 
UR10). Second, we implemented a non-recursive algorithm: we parse the whole cube to find a pixel 
above threshold which is not yet included into a structure ; we tag it and we scan the whole cube several 
times to tag the neighboring pixels of this growing seed which happen to be above threshold, until we 
find no new pixel to attach to this structure ; finally we reiterate to find another pixel not yet included in 
a structure until we fail to find any new pixel above threshold. The second algorithm is much more CPU 
time consuming, but keeps the memory usage constant and is easier to implement. We checked that both 
algorithms identify strictly the same structures for the two implementations in the low resolution cases. 

Like in UR10, all the above quantities are found to exhibit power-law scaling with respect to structure 
linear size Li, with different scaling behavior in the inertial and dissipative ranges. The quantity Xi, 
which could be any of Li, Ri,Vi, Ai, Pi scales as 

Xi oc Lf ■ x 

with different scaling exponents Dx in the inertial and dissipative ranges, while the pdf of Xi scales as 

V(Xi) oc X~ TX 


with different scaling exponents tx in the inertial and dissipative ranges. As an example, the scaling 
relations Pi oc Lf p and V(Pi) oc P~ Tp are shown in figures llR for the structures extracted from Run 12 
(AD - OT), at the peak of dissipation, with a threshold of two standard deviations above mean value. 
The limits of the inertial and dissipation ranges are also shown (as used by UR10, see section 3.4). The 
upper limit of the dissipation range is just below the lower limit of the inertial range, while the lower 
limit of the dissipation range is ~ 2 times the numerical resolution. 


5.3 Comparison with UR10 

In this section we compare the results of the statistical analysis of structures of high dissipation with 
those of UR10. These authors consider pure MHD and study the structures of high Ohmic dissipation 
or high enstrophy 

Suj — ReX 1 ^ 2 , u> = V x u 

The results for the scaling exponents are shown in figure [18] for the case of the OT initial condition. 
The case of the ABC initial condition (not shown) is similar. The exponents are calculated from run 
10 which in terms of initial condition and Reynolds number is similar with run III of UR10. As in the 


25 





Figure 17: Scaling relations Pi oc L^ p from Run 12 (AD - OT), at the peak of dissipation, with a threshold 
of two standard deviations above mean value. The dotted line shows the effect of adopting the slope found 
by Uritsky et ah 2010 instead of our own slope. 
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present paper, the snapshot analyzed is on the peak of dissipation. The structures of high dissipation 
are defined as connected sets of points having values of the Ohmic dissipation two standard deviations 
above the mean value. We keep the same definitions as UR10 for the inertial and dissipative ranges (cf. 
section 3.4). 

Figure 1 18| shows that the agreement between our results and those of UR10 is not completely satisfac¬ 
tory. Although most three-sigma error bars are compatible and the errors are on the same order, there 
remains systematic differences, especially for our pdf exponents which appear to yield shallower pdfs 
than UR10. The least square method used by UR10 to estimate the slope of power-law pdfs is known to 
introduce some bias and the maximum-likelihood estimate method (MLE) should be used instead [see 
Clauset et al. 2009 . We computed the exponents with the MLE, and found them to be very close to 


our least-square values. We turned to explore the effects of some systematics due to the uncertainty on 
the boundaries of the inertial and dissipative range and the size of the bins used to produce the pdfs. 
We varied randomly these parameters with factors in an octave centered on their initial chosen values. 
The excursion of the resulting three-sigma error bars over a thousand of such realizations are plotted as 
thin red error bars on Figure |18| A slight displacement of the inertial or dissipative range boundaries 
incorporates (or leaves out) new data near the edge of the fitting intervals, where their leverage on the 
fitted slope is quite important. The resulting figure shows that such systematics can account for nearly all 
discrepancies with respect to UR10, except for the correlation between the total dissipation of structures 
and their linear size. We illustrate the corresponding discrepancy on Figure [IT] where the dotted red 
line shows the slope followed by UR10 data: this line clearly falls below our data at small scales. This 
suggests that our smallest structure have a higher value of dissipation. This may perhaps be traced back 
to the slightly more refined dealiasing rule which we use. The picture is the same for the ABC runs, 
except for the error bar 0.02 on Da in Run I of UR10 (see their Table II) which is probably a typo as 
the value they quote does not correspond to the scatter displayed in their Figure 3. 

To validate further our results, we compared the computed scaling exponents for the dataset corre¬ 
sponding to the temporal peak of total dissipation of the pure MHD run with ABC or OT initial condition 
(Run 9 and 10) with those computed from a snapshot taken one macroscopic eddy turnover time later, 
in the decay period of the turbulence. In agreement with the results of UR10, we find no statistical 
difference between the two snapshots (one-sigma error bars are compatible). Similarly, we compared the 
exponent values computed based on j 2 with those computed based on w 2 for the dataset corresponding 
to the peak of total dissipation of the same runs (Runs 9 and 10). Again the exponents were seen to be 
compatible within one-sigma error bars, as in UR10. 


5.4 Structures based on total dissipation 

In this section we focus on the statistical analysis of structures extracted based on the total dissipation 
£t = s 0 + £v + £« for both pure MHD and AD MHD. The analysis based on the total dissipation is more 
relevant to the heating of the ISM because all three different types of dissipation can be important heating 
agents. AD has an additional specificity because the ion-neutral drift increases the effective temperature 
of the chemical reactions, but we do not consider this yet in the present work. In the following, we will 
note Dx for the linear size exponents and tx for the probability exponents of a characteristic X. All 
the structures discussed in this section are defined as connected sets of points having values of the total 
dissipation two standard deviations above the mean value. 

In Table [2] we present the results of the structure extraction algorithm. The relative amount of 
dissipation contained by all the detected structures does not depend much on the Reynolds number, but 
the volume filling factor of the structures decreases as the Reynolds number increases. The presence of 
AD results in fewer detected structures than without AD (except for run 13). However, they fill roughly 
the same volume fraction, thus AD structures tend to be larger. This difference between the number 
of different structures in pure MHD and AD increases with the Reynolds number. Figure gives a 
more detailed view of the fraction of total dissipation contained in structures with a value above a given 
threshold as a function of the volume fraction occupied by these structures. The curve rises steeply near 
the origin, so that 30 percent of the dissipation in contained in less than 3 percent of the total volume. 
The steepness at the origin is seen to be mainly due to the Ohmic heating: this is in line with the original 
picture of Brandenburg and Zweibel 1 994] where AD forms sharp features in which Ohmic dissipation 
is favored. 

Figures [20] and [2l] summarize all results on the exponents for the total dissipation two standard 


27 










3.0 

2.5 
2.0 

1.5 
1.0 
0.5 
0.0 

°' 5 D l D r D v D a D p t l r R t v t a t p 

, Dissipative range, OT 


3 

2 

1 

0 


Dl Dp Dy Da Dp ^~r i~ v Ta Tp 



Inertial range, OT 

li ,i *t III 


□ □ Uritsky et al. 
} } Present work 


Figure 18: Scaling exponents for structures extracted based on the Ohmic dissipation (red circles), compar¬ 
ison with results of UR10 (blue squares). Upper panel: Inertial range exponents for our Run 10 (MHD-OT) 
compared with the corresponding Run III of UR10. Lower panel: Dissipative range exponents for the same 
runs. Thick error bars are three-sigma error brackets. Thin red error bars estimate the systematics related 
to the choice of the boundaries of the inertial and dissipative ranges as well as the bin size for the definition 
of the pdfs: see text for details. 
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# 

N 

Re x 

Re a 

7 ^ of structures 

% of volume 

% of total dissipation 

1 

128 

189 

- 

87 

3.90 

23.85 

2 

128 

210 

- 

54 

3.81 

27.29 

3 

128 

217 

100 

72 

3.78 

23.65 

4 

128 

228 

100 

34 

3.90 

26.91 

5 

256 

353 

- 

215 

3.26 

27.32 

6 

256 

377 

- 

233 

3.26 

26.80 

7 

256 

412 

100 

144 

2.86 

26.40 

8 

256 

444 

100 

153 

2.98 

26.53 

9 

512 

591 

- 

790 

2.70 

31.30 

10 

512 

604 

- 

1166 

2.77 

30.33 

11 

512 

756 

100 

375 

2.17 

29.21 

12 

512 

750 

100 

418 

2.57 

29.48 

13 

512 

640 

10 

1167 

2.94 

28.01 

14 

512 

927 

10 

378 

1.75 

22.00 


Table 2: Results of the structure extraction algorithm for all runs and structures defined as connected sets 
of points having a value of the total dissipation two standard deviations above mean value. 



Figure 19: We look at the subset of pixels above a given threshold of total dissipation in run 12 (AD-OT, 
Re a = 100) at the dissipation peak. For each value of the threshold, we plot the fraction of the total energy 
dissipation on this subset versus the volume of this subset (black curve). We also give the fraction of the 
total dissipation on this subset for each nature of dissipation (red: Ohmic, green: viscous, blue: AD). 
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Figure 20: Comparison of scaling exponents with three-sigma error bars between pure MHD (red circles) 
and AD MHD (green squares Re a = 100 and blue triangles Re a = 10) - ABC Runs 9,11,13 
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Figure 21: Same as Figure |20| but for the OT Runs 10,12,14. 
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Figure 22: Scatter plots of the ratio e 0 /e t versus the ratio e a /st for AD-OT runs 12 and 14. 


deviations above the mean. Although most three-sigma error bars are compatible, in particular for the 
Dx exponents which are unchanged with AD, some systematic differences exist for the pdfs exponents tx- 
The pdfs exponents are in general steeper in the dissipation range: AD seems to favor more fragmented 
structures in the dissipative range. This confirms the tendency for more intermittency with AD that 
was suggested by the structure functions analysis. However, we see no clear cut tendency in the inertial 
range. If we discard the strong AD results (. Re a = 10), the inertial range shows a behavior opposite 
from the dissipative range (shallower pdfs slope, ie: larger structures are favored when AD is present). 
However if we now look only at the MHD runs and the Re a = 10 runs, the inertial range sees no change 
in the tx exponents for the OT runs, but bigger exponents for the ABC runs, in agreement with the 
dissipative range and contrary to the Re a = 100 runs... 

This complicated picture might perhaps not be genuine, as the huge systematics experienced for the 
comparison with UR10 show. However, the dependence of the intermittency statistics on the initial 
conditions (seen both in the structure pdfs slopes and in the structure functions exponents) points at 
their difference in magnetic helicity content. In the OT initial condition, the initial value of magnetic 
helicity is almost zero, and the equation of magnetic helicity evolution is unchanged by the inclusion of 
the AD term. As far as the effects of viscosity and resistivity are neglected, both the pure MHD and AD 
MHD solution evolve under the same constraint of zero magnetic helicity. A zero value for the magnetic 
helicity is an important constraint because it implies statistical reflection invariance, a property which 
the ABC runs will not share. In the ABC initial condition, the constraint of very low cross-helicity is 
broken by AD, which provides a source term in the cross-helicity equation. 

The above analysis does not give any information on the relative amount of Ohmic, viscous and AD 
dissipation contained within each structure. To answer this question, we show in figure [22] scatter plots of 
the ratio of total Ohmic dissipation to total dissipation £ 0 /st versus the ratio of total AD dissipation to 
total dissipation e a /st in the two OT runs 12 and 14. In all AD cases there is a tendency for the structures 
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to cluster close to the line of zero viscous dissipation (dashed line in figure 22). This shows that within 
the structures of high total dissipation, viscous dissipation is relatively less important. The relative 
values of Ohmic and AD dissipation span however the whole spectrum, in contrast to the impression 
given by our RGB slices (see figures [7] to [9| : this could be a genuine difference between the extreme 
dissipation events and the bulk of the dissipation shown on the RGB figures, or the extraction algorithm 
of connected structures could merge nearby sheets with different dissipation natures. We checked on a 
pixel by pixel scatter plot similar to Figure [22] that it is indeed a genuine difference. For the strongest 
AD runs, though, the intense dissipation structures tend to be predominantly due to AD heating (top 
left corner in Figure 22). 


6 Concluding remarks 


We performed a 3D numerical study of the structures of high dissipation in MHD turbulence, with the 
inclusion of ambipolar diffusion. At the Reynolds numbers studied, the total dissipation due to viscosity, 
resistivity and AD are of comparable magnitude. 

Kinetic and magnetic energy spectra show that ambipolar diffusion enhances the turbulent energy to 
small scales at the expense of intermediate scales. This agrees with the idea of |Brandenburg and Zweibel] 
[1994] that AD can sharpen the magnetic gradients, but the effect is not strong enough to increase the 
total Ohmic dissipation rate in our simulations. Previous authors Li et al. [2008] and Oishi and Mac Low 


[2006] have examined the case of driven two-fluids compressible turbulence with a mean magnetic field 


Oishi and Mac Low 2006 find no effect on the slope of the spectrum (for 
Li et al. 2008 find that AD steepens it, although they don’t display the 


and they find different results: 

Re a = 2.5 and Re a — 5) while 
results for their runs with Re a = 12 and Re a = 120, the only ones which have £ a in the computed range 
of scales as in our study. It should be noted that both these papers neglect Ohmic diffusion, and rely 
on truncation errors only to reconnect the field: the present work is the first to account for AD in the 
presence of controlled Ohmic and viscous dissipation. This is important because the dissipation physics 
can be quite different from the numerical dissipation as was demonstrated by |Fromang and Papalo izou| 
[2007] , Fromang et al. [2007] in magneto-rotational turbulence. 

As in Oishi and Mac Low 2006 , we fail to detect a significant change of regime in the spectra at the 
expected AD length scale £ a , but in our simulations it happens at a greater length scale. This length 
scale £ a is predicted from the balance between the moduli of the Fourier coefficients for the inertial e.m.f. 
and the AD e.m.f.. We would underestimate £ a if AD was more coherent in time than the advection of 
the held, and so we conjecture that AD terms have a greater coherence time than the inertial terms. 

In our simulations, we observe that AD shuts off the Lorentz force at small scales: this shifts the 
peak of the AD heating power spectrum to larger scales. The position of that peak defines a scale which 
seems to match the characteristic thickness of the sheets where AD heating is strong. This scale £* might 
be revealed by the characteristic chemistry of AD heating where neutral-ion endothermic reactions are 
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favored (e.g.: CH + or SH + formation 
of various environments of the ISM. We identify the integral length scale and the r.m.s 


2009). Table [ 3 ] sums up the characteristics 

velocity in 

our simulations to their corresponding physical values in each considered ISM components to apply our 
results. Although we find that Re” 1 in the ISM varies from 10 -3 to 10 -2 , we use £* s= 0.6 as measured 
in our AD simulations with Re~ x — 0.01 to estimate the scale lo£* of the AD diffusion heating. The AD 
dissipation heating is always about a fourth of the typical scale whereas lo.£ a can be much smaller. 

The qualitative picture of the dissipation field suggests that it is dominated by intermittent sheet-like 
structures which alternate with large voids of low dissipation. The sheets of various dissipative nature 
(Ohmic, viscous or AD) appear to be clearly separated, except for the highest dissipation rates, where 
viscosity fades out and Ohmic and AD heating blend. AD heating sheets are often seen to sandwich much 
thinner regions of strong Ohmic or viscous dissipation, as in the simple case studied by |Brandenburg| 
and Zweibel 1994 . The high degree of intermittency is confirmed by the computation of the structure 
function exponents for the velocity and the magnetic field, as well as by the pdf of the total dissipation, 
which exhibits a log-normal core and a strong power-law tail for high values of the dissipation. We 
compared the statistics of the structures of strong dissipation with those of UR10 and we obtain only 
marginal agreement, probably because of the systematics linked with the definition of the inertial and 
dissipative ranges. The statistical analysis of structures of high total dissipation reveals the highly 
intermittent nature of the dissipation held, as more than 30% of the dissipation takes place in less than 
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Table 3: Characteristics of various components of the ISM. Dimensions are recovered from our simulations 
by assuming L ~ 2.5, ~ 0.6 and uo ~ V?>.U a where U a is the line-of-sight r.m.s. Alfven velocity. The 

quantity 7 pi is computed by assuming the ions are essentially C + ions with a number density rii = 10 - 4 nn. 




CNM 

molecular clouds 

low-mass dense cores 

Density 

nn(cm 3 ) 

30 

200 

10 4 

Length scale 

L.l 0 (pc) 

10 

3 

0.1 

r.m.s. velocity 

U.uo/\/3 (km/s) 

3.5 

1 

0.1 

Alfven velocity 

U a = u 0 /\/3 (km/s) 

3.4 

2 

1 

AD Reynolds number 

-Re" 1 = I /7 Pi/(to) 

1.2 10“ 2 

3.6 10“ 3 

1.1 10“ 3 

AD heating length 

t a .l 0 (pc) 

2.4 

0.72 

0.024 

AD dynamical length 

£ a .l 0 = 2nl 0 Re- 1 .U 2 a /U 2 (pc) 

0.28 

0.10 

0.026 


3% of the volume. No significant difference in the scaling laws between pure MHD and AD MHD was 
found, but the slope of the power-law pdfs is affected in the dissipative range, with a statistical preference 
towards more fragmented structures. 

In future work, we intend to make progress towards a more realistic picture of the ISM, relaxing the 
incompressible hypothesis, with the final aim of including realistic cooling. We also hope our statistical 
results will provide new ways to approach the observed characteristics of the ISM, intermediate between 
direct post-processing of 3D numerical simulations and the building of line of sights with elementary 
models such as shocks or vortices [Godar d et al. 2009] . 
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